https://github.com/prioriactions/teaching

Reproducibility, Academic use, …

Source: https://medium.com

Three questions: Planning objective, Optimization objective, and Intensity of threats.

Recovery vs Conservation

Binary vs Continuous

Bolch, E.A. et al. (2020). Remote Detection of Invasive Alien Species

Binary vs Continuous

Source: https://ib.bioninja.com.au

Data that defines the spatial prioritization problem (planning units data, feature data, threats data, and their spatial distributions).

5 + 2 (optional)

inputData(
  pu,
  features,
  dist_features,
  threats,
  dist_threats,
  sensitivity = NULL,
  boundary = NULL
)

What if I have my data in another format?
A: There are functions that allow you to import and convert data to data.frame format

Examples:

##from .txt to data.frame
df <- read.table('myfile.txt',sep='\t')

##from .csv to data.frame
df <- read.csv("myfile.csv")

Source: https://libguides.umn.edu

data.frame

Specifies the planning units (PU) of the corresponding instance and their corresponding monitoring cost and status. Each row corresponds to a different planning unit. This file is inherited from the pu.dat in Marxan.

## id monitoring_cost status
## 1                2      0
## 2                2      0
## 3                2      0
## 4                2      0
## 5                2      0
## 6                2      0

id

Game, E. T. y H. S. Grantham. (2008). Manual del Usuario de Marxan: Para la versión Marxan 1.8.10

monitoring_cost

Source: https://marxansolutions.org

status (optional)

• free (0)
• locked-in (2)
• locked-out (3)

Source: https://marxansolutions.org

data.frame

Specifies the conservation features to consider in the optimization problem. Each row corresponds to a different feature. This file is inherited from the spec.dat in Marxan.

##id target_recovery     name
## 1              11 feature1
## 2              16 feature2
## 3               8 feature3
## 4               9 feature4

id (Species)

Source: https://gefespeciesamenazadas.mma.gob.cl

id (Ecosystem Services)

Source: https://ecology.fnal.gov

id (Habitats)

Source: https://www.asemafor.cl/areas-protegidas-en-chile

recovery_target

conservation_target (optional)

data.frame

Specifies the spatial distribution of conservation features across planning units. Each row corresponds to a combination of planning unit and feature.

#pu feature amount
#1        3      1
#2        3      1
#3        3      1
#4        3      1
#5        3      1
#6        3      1

amount (binaries)

Source: https://calumma.typepad.com

amount (continuous)

Source: https://ib.bioninja.com.au

data.frame

Specifies the threats to consider in the optimization exercise. Each row corresponds to a different threats.

#>id    name blm_actions
#> 1 threat1           0
#> 2 threat2           0

blm_actions

Source: https://docs.marxanweb.org

data.frame

specifies the spatial distribution of threats across planning units. Each row corresponds to a combination of planning unit and threat.

#>pu threat amount action_cost status
#> 8      2      1           2      0
#> 9      2      1           2      0
#>10      2      1           2      0
#>11      1      1           3      0
#>11      2      1           4      0
#>12      1      1           3      0

data.frame

specifies the spatial distribution of threats across planning units. Each row corresponds to a combination of planning unit and threat.

• action_cost
  numeric cost of an action to abate the threat in each planning unit.

• status (optional)
  integer value that indicate if each planning unit should be available to be selected (0), locked-in (2) as part of the solution, or locked-out (3) and excluded from the solution. :::

amount (continuous)

Source: https://ib.bioninja.com.au

Optional

data.frame

specifies the sensitivity of each feature to each threat. Each row corresponds to a combination of feature and threat.

#>feature threat
#>      1      1
#>      2      1
#>      3      1
#>      4      1
#>      1      2
#>      2      2

Optional

data.frame

specifies the sensitivity of each feature to each threat. Each row corresponds to a combination of feature and threat.

• c (optional)
  numeric minimum probability of persistence of a features when a threat reaches its maximum intensity value.

• d (optional)
  numeric maximum probability of persistence of a features in absence threats. :::

a, b, c and d

a, b, c and d

Optional

data.frame

Specifies the spatial relationship between pair of planning units. Each row corresponds to a combination of planning unit.

#>id1 id2 boundary
#>  1   1        0
#>  2   1        1
#>  3   1        2
#>  4   1        3
#>  5   1        4
#>  6   1        5

Minimize Costs (minimizeCosts) vs Maximize Benefits (maximizeBenefits)

\(x_{ik}\) is the decisions variable that specifies whether an action to abate threat \(k\) in planning unit \(i\) has been selected \((1)\) or not \((0)\).

\(c_{ik}\) is the cost of the action to abate the threat \(k\) in the planning unit \(i\).

\(b_{is}(x)\) is the benefit of the feature \(s\) in the planning unit \(i\) (ranging between \(0\) and \(1\)).

\(t_s\) is the recovery target for feature \(s\).

\(f(x)\) is the function of connectivity penalty.

Minimize Costs (minimizeCosts) vs Maximize Benefits (maximizeBenefits)

\(x_{ik}\) is the decisions variable that specifies whether an action to abate threat \(k\) in planning unit \(i\) has been selected \((1)\) or not \((0)\).

\(c_{ik}\) is the cost of the action to abate the threat \(k\) in the planning unit \(i\).

\(b_{is}(x)\) is the benefit of the feature \(s\) in the planning unit \(i\) (ranging between \(0\) and \(1\)).

\(B\) is the budget available.

\(f(x)\) is the function of connectivity penalty.

Create an optimization model for the multi-action conservation planning problem, following the mathematical formulations used in Salgado-Rojas et al. (2020).

problem(
  x,
  model_type = "minimizeCosts",
  budget = 0,
  blm = 0
)

Solves the optimization model associated with the multi-action conservation planning problem. This function is used to solve the mathematical model created by the problem() function.

solve(
  a,
  solver = "",
  gap_limit = 0,
  time_limit = .Machine$integer.max,
  solution_limit = FALSE,
  cores = 2,
  verbose = TRUE,
  name_output_file = "output",
  output_file = TRUE
)

solver (gurobi, cplex, symphony)

gap limit (percentage)

Source: https://www.math.uwaterloo.ca

time limit (seconds)

solution_limit (First solution?)

verbose (Information on screen?)

Returns the spatial deployment of the actions for each planning unit of the corresponding solution.

getActions(x, format = "wide")

Output example,

##   solution_name pu 1 2 conservation connectivity
## 1           sol  1 0 0            0            0
## 2           sol  2 0 0            0            0
## 3           sol  3 0 0            0            0
## 4           sol  4 0 0            0            0
## 5           sol  5 0 0            0            0
## 6           sol  6 0 0            0            0

Returns the total benefit induced by the corresponding solution. The total benefit is computed as the sum of the benefits obtained, for all features, across all the units in the planning area.

getSolutionBenefit(x, type = "total")

Output example,

##   solution_name feature benefit.conservation benefit.recovery benefit.total
## 1           sol       1                    0               11            11
## 2           sol       2                    0               16            16
## 3           sol       3                    0               10            10
## 4           sol       4                    0                9             9

Provides the sum of costs to actions and monitoring applied in a solution.

getCost(x)

Output example,

##   solution_name monitoring threat_1 threat_2
## 1           sol         61       20       65

Package prioriactions can be found at CRAN, so it can be installed using:

install.packages("prioriactions")

Latest stable versions can be downloaded and installed from GitHub as follows (package remotes should be installed first):

if (!require(remotes)) install.packages("remotes")
remotes::install_github("prioriactions/prioriactions")

So, we load the prioriactions package.

# load package
library(prioriactions)

(Available in the Get Started vignette on the prioriactions website)

• This example contains 100 planning units, 4 features and 2 threats.

• The distribution of features and threats can be plotted on a grid of 10 x 10 units.

• prioriactions contains this simulated example inside the setup files. You can extract it by the data() function.

The monitoring cost values ranging from 1 to 10 and all status of 0 (not locked).

# load planning unit data from prioriactions
data(sim_pu_data) #To load simulated data

#head(sim_pu_data)

A RasterLayer object can be used to present this spatial information. The pixel values correspond to the monitoring costs of each planning unit.

library(raster) #To plot rasters

r <- raster(ncol=10, nrow=10, xmn=0, xmx=10, ymn=0, ymx=10)
values(r) <- sim_pu_data$monitoring_cost
plot(r)

Contains information about the features such as its id and targets (mandatory when minimizing costs).

# load features data from prioriactions
data(sim_features_data)

#head(sim_features_data)

Contains information on the spatial distribution of these features across planning units.

# load features data from prioriactions
data(sim_dist_features_data)

#head(sim_features_data)

To plot the spatial distribution of the first feature,

# load amount of features data
features <- reshape2::dcast(sim_dist_features_data, 
                            pu~feature,
                            value.var = "amount", 
                            fill = 0)

values(r) <- features$`1`
plot(r)

Provides information about the threats such as their id and name.

# load threats data from prioriactions
data(sim_threats_data)

Provides information on the spatial distribution of these threats

# load threats data from prioriactions
data(sim_dist_threats_data)

Indicates which features is sensitive to what threat.

# load threats data from prioriactions
data(sim_sensitivity_data)

Provides information on the spatial relationship between PU and they are presented in long format.

# load boundary data from prioriactions
data(sim_boundary_data)

After having loaded our data, we will now create the data object through the inputData() function.

# create conservation problem
b <- inputData(pu = sim_pu_data, 
               features = sim_features_data, 
               dist_features = sim_dist_features_data, 
               threats = sim_threats_data, 
               dist_threats = sim_dist_threats_data, 
               sensitivity = sim_sensitivity_data, 
               boundary = sim_boundary_data)

# print problem
print(b)

## Data
##   planning units: data.frame (100 units)
##   monitoring costs:     min: 1, max: 10
##   features:       feature1, feature2, feature3, feature4 (4 features)
##   threats:        threat1, threat2 (2 threats)
##   action costs:   min: 1, max: 10

It is advisable to determine the maximum benefit achievable for each conservation feature through the getPotentialBenefit() function

# get benefit information
getPotentialBenefit(b)

# create optimization problem
c <- problem(b, model_type = "minimizeCosts")

## Warning: The blm argument was set to 0, so the boundary data has no effect

## Warning: Some blm_actions argument were set to 0, so the boundary data has no
## effect for these cases

# print model
print(c)

## Optimization Problem
##   model sense: minimization
##   dimensions:  284, 396, 12.656 kB (nrow, ncol, size)
##   variables:   396

To solve our model, we need an optimizer. In this case, we will use the Rsymphony library.

# solve optimization problem
d <- solve(c, solver = "symphony", verbose = TRUE, output_file = FALSE, cores = 2)

# print solution
print(d)

## Solution overview
##   name: sol
##   objective value: 738
##   gap:  0%
##   status:  Optimal solution (according to gap tolerance: 0)
##   runtime: 0.017 sec

getActions()

# get actions from solution
actions <- getActions(d, format = "wide")

# print first six rows of data
# head(actions)

getActions()

values(r) <- actions$`1`
plot(r)

We go back to step 2 and repeat.

# create optimization problem
c2 <- problem(b, model_type = "minimizeCosts", blm = 10)

## Warning: Some blm_actions argument were set to 0, so the boundary data has no
## effect for these cases

# print model
print(c2)

## Optimization Problem
##   model sense: minimization
##   dimensions:  29984, 10296, 883.856 kB (nrow, ncol, size)
##   variables:   10296

# solve optimization problem
d2 <- solve(c2, solver = "symphony", verbose = TRUE, output_file = FALSE, cores = 2)

# print solution
print(d2)

## Solution overview
##   name: sol
##   objective value: 863
##   gap:  0%
##   status:  Optimal solution (according to gap tolerance: 0)
##   runtime: 22.847 sec

# get actions from solution
actions2 <- getActions(d2, format = "wide")

values(r) <- actions2$`1`
plot(r)

• We divided the whole catchment (71,630 km2) into 2316 sites (i.e., sub-catchments)

• We sourced the distribution of 45 fish species in the Mitchell river catchment as our conservation features.

• We considered four major threats to freshwater fish species in the catchment: water buffalo (Bubalis bubalis), cane toad (Bufo marinus), river flow alteration (caused by infrastructure for water extractions and levee banks) and grazing land use.

https://prioriactions.shinyapps.io/prioriactions/